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Abstract For most hierarchical triple stars, the classical double two-body model of 
zeroth-order cannot describe the motions of the components under the current observa¬ 
tional accuracy. In this paper, Marchal’s first-order analytical solution is implemented and 
a more efficient simplified version is applied to real triple stars. The results show that, for 
most triple stars, the proposed first-order model is preferable to the zeroth-order model 
both in fitting observational data and in predicting component positions. 
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1 INTRODUCTION 

A hierarchical triple star is composed of a close binary and a distant third component. About one 
thousand stars of this kind are contained in the latest on-line version of The Multiple Star Catalog 
jTokovinin 19971) . In these systems, the primary components are usually bright. Bright stars are useful 
in many aspects (e.g. lUrban & Seidelmann 20141) . Though a set of isotropic and dense stars is crucial 
for some applications such as navigation, the stars with nearby companions are usually excluded. This 
is the case for the Hipparcos Celestial Reference Frame, as recommended in lAU resolution B1 (2000fl 
For triple stars, the problem lies mainly in that the primary positions generally cannot be predicted 
accurately by the almost exclusively used model, namely the classical double two-body model. 

Hierarchical triple stars are also of great interest in stellar physics and galactic astronomy, due 
to the fact that their dynamical evolution is important to both stellar and galactic evolutions (e.g. 
Binney & Merrifield 1998| IValtonen & Karttunen 20061 lAarseth 20031) . Moreover, these systems are 
often studied in terms of stability of the general three-body problem (e.g. IMarchal & Bozis 19821 
|Li, Fu & Sun 2009| ). In some case studies, the results are sensitive to the mass parameters and the ini¬ 
tial conditions (e.g. lOrlov & Zhuchkov 2005l l. the accuracies of which are limited again by the double 
two-body model used in fitting observations (e.g. ILiu et al. 20091 1 . 

As a zeroth-order solution of the hierarchical three-body problem, the double two-body model has 
the advantage of being analytical and simple. The existing first-order analytical solutions are more ac¬ 
curate. The former one is still dominantly used, while the latter ones, as far as we know, remain little 
used in fitting observations. In this paper, the first-order solution by Marchal is efficiently implemented. 
This is achieved mainly by making some simplified modifications and high order approximations to 
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Marchal’s solution. In the context of fitting observations of triple stars, we call Marchal’s solution and 
the double two-body solution, respectively, the M-model and the K-model. 

In section 2, the M-model is implemented. In section 3, the improvement in accuracy of M-model to 
K-model is statistically discussed with a set of sampling triple stars. In section 4, a simplified M-model 
is given and applied to real triple stars. Concluding remarks are given in the last section. 

2 AN IMPLEMENTATION OF M-MODEL 

Consider a hierarchical three-body problem in an inertial coordinate system {0—xyz}, where O is the 
center of mass and the z-axis parallel to the total angular momentum C. Denoting the masses of the 
inner two bodies by mi and m 2 , and the mass of the third body by m 3 , we will use the following 
mass-dependent parameters. 


mt = mi + m 2 + m 3 , 


rrii 


mim 2 


mi+m 2 ’ 


A = 


1721+7712 ’ 


Q _ G ( 1711 + 7712 ) 

“ mt 


_ {mi+m2)m3 
° mt ’ 

o _ G^(mi+m2)^rra3 

(mim 2 mt)^ 


where G is the gravitational constant. Let r be the position vector of m 2 relative to mi, and R the 
position vector of m 3 relative to the center of mass of the binary. The ratio e = -^ = is a small 
quantity. 

The Delaunay variables as expressed in terms of the ordinary orbital elements (a, e, i, ui, il, M) are 


Gi — m^ G(^m,i -t- m2)u2, Qi — C*\/l , R-i — Qt cos , 

ii = Afj, gi = uji , hi = flj, 

Gq — GiTitCio , Qo — 1 Cq , Ifo — Qo cos %Q 

to — tVio , go — iOo , ho = flo , 


where the subscripts i and o indicate the inner and outer orbits, respectively In these variables, the 
Hamiltonian up to the first order in ^ formally written as 

H = H{£i,gi,£o,GoJi,gi,to,go,Ri + Ro,ho- K) 


~ Hoi + Hqo + Hi 

^ _^ ^ + (1 - e^cosL;,)^ t _ 2\ f 

2£2 2£l 2£lil-eo cos Eof^ ^\CoJ ’ 


( 1 ) 


where $ = ^{£i,Gi,£o,GoJi,giJo,go,'Hi + Ra.K - hi) = ^,Ei = Et{£i,Gi;ti) and K = 
Eo{£o, Go', to) are the eccentric anomalies of the inner and outer orbits, respectively. 

In eq.([T]i, TLo + Hi and ho — hi are understood as two single canonical variables conjugating re¬ 
spectively to the negligible hi and TLo- And so, they are constants that can be calculated from the initial 
conditions. The standard way to calculate the two negligible variables is by quadrature, after all the other 
degrees of freedom are integrated. But in the present context, we have as consequences of the integral 
of angular momentum 


nt + no = c = \ci ho-ht = TT, no = 


C^ + G^o-Gf 
2C 


Therefore, only hi needs to be calculated by quadrature. Because of the short-period terms in the inte¬ 
grand, the numerical quadrature is time-consuming. It is then preferable not to follow the standard way 
and decouple only {TLo, ho — hi) from the other degrees of freedom at this stage. 

For the system defined by the Hamiltonian eq.([T]) with ho — hi = tt, a first-order integrable system 
can be achieved by the Von Zeipel transformation (e.g.|Harrington 19^ 1969. IMarchal 197^ 1990). In 
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the resulting canonical variables {Ci,Qi,Co,Go,C]£i,gi,£oTgo,hi), called long-period Delaunay 
variables, the new Hamiltonian can be written as 


where 


H = H{Li,Qi,Lo,Qo,C,gi) 

= Hqi -I- Hqo + Hi 

A Po , Pii'iz - £>)Co ( 

£^o 


2L] 


2Ll 






2 - 


- g| - Ql 

2 .GIG 0 




1 - 


C" - gf - Gl 
2 .G 1 G 0 


sin^ gi. 


( 2 ) 


(3) 


In this time-independent Hamiltonian of five degrees of freedom, there are four negligible variables 
£i,£o,go,hi. Their conjugate variables £/, £ 0 , Go and C, together with the total energy H and 
z = z{H,£[,Co,Go) as given by solving eq.(|2|i, are constants known from initial conditions. This 
confirms the integrability of the transformed Hamiltonian system. 

The differential equations for Gi and gj, the variables of the only non-negligible degree of free¬ 
dom, can be integrated simultaneously. But to be more efficient, we first integrate the equation for Gi, 

q2 

decoupled from gj by using eq.Q. In terms of € (0,1), this equation writes 

L.J 

i = (4) 

where, with A = and B = 

Pi{x) = B'^x^ — 2(1 + AB)x + z + A^, 

P 2 (x) = 4B^x^ - (5B^ + SAB + 3)x^ + (4^^ + lOAB - z + 5)x - 5A^. 


From the necessary condition Pi{x)P 2 {x) > 0, Marchal (1990) pointed out that x oscillates between 
two neighbouring roots, Xa G (0,1) and Xb G {xa, 1), of Pi{x)P 2 {x). To be specific, the function x{t) 
defined in eq.® changes its sign from negative to positive at Xa, and the opposite is true at Xb- 

The difficulty in integrating eq.® caused by this unfavorable feature of the right-hand side can be 
avoided. For this, we introduce a continuously changing angular variable 9, for which mod(27r) is not 
allowed, by the following variable substitution x = Xa + (xb — x a) sin^ 9. 

Let 173 , fT 4 , (T 5 be the other three roots of Pi {x)P 2 (x). We have 


where 


(0)= _1_ 

d0 ^ y/'l — ci sm^(e)+C2 sin'‘(0) —C3 sin®( 6 l) ’ 


3 /3i£j 




Ba-t, a = > 0 : 


di + c?2 + da, C2 — did2 + dids -I- d2d3, 
Xb - Xa 


Cl 

C 3 = did 2 d 3 > 0 , 


dj — 


^J+2 


-,(j = 1,2,3). 


(5) 


Given the initial condition {to, do), the value of 9 at any time t can be obtained from an iterative method. 
And, given 9, Gi{> 0) can be calculated from the defining formulae of 9 and x. 

As I sinp/(f)| can be solved from eq.@, the key to determining gj is its quadrant. Let n be the 
biggest integer no greater than 29/tt. The quadrant of gi{t) can be deduced from the type of motion, 
gi (0) and 9. Depending on the initial conditions, there are three types of motion. 
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Type 1: P 2 {xa) = 0 and P 2 {xb) = 0. In this type of motion, gi oscillates around ^ or periodi¬ 
cally. In the case of sin((;7(0)) > 0, gi(t) is in the first quadrant if n is odd and the second quadrant if 
n is even. In the other case, gi{t) is in the third quadrant if n is odd and the fourth quadrant if n is even. 

Type 2; P 2 {xa) = 0 and Pi{xb) = 0. In this case, gj always increases as time grows. The gi(t) is 
in the same quadrant as [0„, + ^), where On = if g/(0) is in the same quadrant as [—f, §) , 

and On = if gi{0) is in the same quadrant as [f, %-) . 

Type 3: P\{xa) = 0 and P 2 {xb) = 0. The gj always decreases as time goes by. The gi{t) is in the 
same quadrant as (0„ — where On = (l — ^) tt if ^/(O) is in the same quadrant as (0 , tt], and 

if 5/(0) is in the same quadrant as (— tt, 0 ]. 

The other four angular variables can be obtained by quadrature, 

£i{t) =£7(0) + |if + /gVi(x(t^))Ii(t?)dt?, £o(f) =fo(0) +Iff+1^(5-3z)f, 

5o(t) =5o(0)+/^„F2(x(t9))Ii(^)dt^, hi{t) = hi{Q)+jlF^{x{m^m^, 

where 


Fi{x) 

F2{x) 

Fsix) 


1 

B^cr 

5-32 

2B(T 


('r _ 5^ I x{z-2) + {A-Bx)^ 

3''“^ 2{l-x) 

+sik 'r-TlYf;* [1 + - B*)i, 

C (z-x){A-Bx) 


' 2B^cr So x-{A-Bxy 


If the first-order long-period solution is gotten, one can make inverse transformations of the solution 
to the original coordinate system. 


3 COMPARISON BETWEEN M-MODEL AND K-MODEL 

In order to compare the accuracy of different models in calculating the observational quantities, it is 
necessary to do a numerical experiment. For the time being, we are interested in only the systems with 
negligible 2nd-order perturbations. Therefore we generated 1000 systems, which satisfy |iT 2 |/|ffoi + 
F[[)o+Hi\ < O.Olin [—100, lOOjyears, andiJ 2 is the second-order perturbation term in the Hamiltonian 
©. This time span is used because the practical cycle of a star catalog is usually less than one hundred 
years. As expected, for some of the generated systems, especially for the systems with large periods 
and high eccentricities of the outer orbits, the first-order averaged perturbations are too large. For such 
a case, M-model fails to be the first-order model. We just consider the samples that satisfy 

|iJi/iJoz| < 0.1, |iJi/iJoo| <0.5, (6) 

during [—Pt, Pt] years, where Pt > max( 100 , Pq), and Pg represents the initial period of the outer orbit. 
Nearly 90 samples are excluded by condition eq.®. In addition, Delaunay elements are not effective in 
describing the orbits that are near circular, near parabolic or near the reference plane, and M-model is 
not suitable to be used in coplanar motion. If there are very small divisors, the implicit Zeipel transfor¬ 
mations can not be solved by the iterative method. Another ^ 40 samples are excluded, and 870 samples 
remain. The remnant samples are used to do a numerical experiment to check the accuracy of M-model 
compared with K-model. 

We calculate the positions of three bodies in the center-of-mass frame during the [—100,100] years 
by M-model and K-model, respectively. As a comparison standard, these positions are also calcu¬ 
lated by the numerical solution {N-model for short). Denote the root-mean-squared errors (RAISE) of 
the 9-dimensional vectors of M-model relative to those of N-model by dM, and the RMSE of the 9- 
dimensional vectors of K-model relative to those of N-model by dx- When (r/R)^ <C {mi +m 2 )/mt, 
generally dM/dx 1, as is shown in Fig. [T] Fig. [T] shows that M-model is apparently better than 
K-model in accuracy when the abscissa is smaller than —1.4. When the abscissa is greater than —1.4, 
Fig. [Ureveals that for most samples the M-model is still more accurate than K-model. 
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Fig. 1 The abscissas on the x-axis are calculated in [—Pt , Pt] years. The abscissa of the dashed 
line is —1.4. Circular points represent the samples that satisfy dm/dK < 0.9, while square 
points represent the samples that satisfy djvr /dfc > 0.9. There are 798 circular points and 72 
square points. 

For a few samples which are at the up-right quarter of Fig. [T] the accuracy of M-model is not as 
good as that of K-model. The phenomena can be explained by the perturbations and the improper use of 
the Delaunay elements. 

There is one sample whose ordinate is apparently greater than 0.5 in Fig.[T] We found that the outer 
orbit of this sample has a very large period and high eccentric. The max (^rni+rn 2 )/mt re^Uy small 
during the considered [—100,100] years, and K-model is very approximate to N-model. While M-model 
considers the averaged perturbations which are much greater. We calculated max \ H — H^i — Hqo — Hi \ 
in [—Pt,Pt] years and max \H — Hm — Hqo \ in [—100,100] years. The former is more than 1000 times 
of the latter, and this supports that M-model is not a first-order model in such cases. 

As the abscissas of samples represented by squared points are not sufficiently small (bigger than 
— 1.4), the inaccuracies caused by small divisors cannot be ignored. For some samples represented by 
squared points in Fig. [1] the detailed reasons are complex and uncertain currently. In all, M-model is 
better than K-model in accuracy for ~ 80% of the samples, and can be credibly applied when the 
abscissa is smaller than —1.4. 

4 THE APPLICATION 

Simplifications of M-model can be made according to the results of the numerical experiment. In eq.(|5]l, 
x{9{t)) can be solved efficiently by an approximation. Generally Ii (r^) can be written 


Xi{9)=X2{d) + [Xi{d)-X2{d)], 


( 7 ) 


where X 2 (d) can be defined as 



if cf - 4 c 2 > 0, 1 - Cl + C2 > C3 > 0, C2 > 0, 
if Cl - 4 c2 < 0,1 - Cl > |c 2 - C 3 I > 0 , 
















6 


X.B. Xu, F. Xia & Y.N. Fu 


Table 1 The application results of the 25 observed triple stars during the time span from 
1900.0 to 2100.0. 


system name perturbation order 

(WDS) (log,„ 

dM 

ii) (AU) 

dMC 

(AU) 

dK 

(AU) 

Improvement 

(l°gio 

Type 

(1/2) 

00325-1-6714 

-1.52 

2.77E-2 

2.78E-2 

0.47 

-1.23 

2 

01148-1-6056 

-4.64 

6.11E-7 

4.96E-4 

1.96E-3 

-0.59 

1 

02022-1-3643 

-1.56 

0.013 

0.013 

0.23 

-1.25 

1 

03082-1-4057 

-4.36 

7.05E-4 

1.43E-2 

8.32E-2 

-0.76 

2 

04142-1-2812 

-4.13 

4.78E-5 

l.lOE-4 

0.10 

-2.96 

1 

04400-1-5328 

-1.53 

0.119 

0.119 

0.96 

-0.91 

2 

06262-1-1845 

-7.69 

2.13E-7 

2.63E-6 

6.01E-5 

-1.36 

2 

07201-1-2159 

-7.35 

7.71E-9 

7.58E-7 

1.23E-5 

-1.21 

2 

10373-4814 

-2.88 

2.60E-4 

2.41E-3 

2.72E-2 

-1.05 

2 

10373-4814 

-2.77 

3.49E-4 

4.73E-3 

3.60E-2 

-0.88 

2 

11308-1-4117 

-6.22 

1.23E-7 

4.96E-6 

4.06E-4 

-1.91 

2 

12108-1-3953 

-1.64 

0.180 

0.180 

0.99 

-0.74 

2 

12199-0040 

-3.24 

1.31E-3 

3.26E-3 

0.18 

-1.74 

2 

15183-1-2650 

-1.76 

0.014 

0.014 

0.12 

-0.93 

2 

16578-1-4722 

-2.39 

6.97E-4 

8.63E-4 

1.66E-3 

-0.28 

2 

17539-3445 

-7.14 

4.58E-7 

2.47E-5 

9.87E-5 

-0.60 

2 

19155-2515 

-4.08 

2.06E-5 

2.03E-4 

1.89E-2 

-1.97 

1 

20396-1-0458 

-1.45 

7.17E-2 

7.17E-2 

1.30 

-1.26 

1 

20475-1-3629 

-2.15 

1.19E-3 

1.19E-3 

5.27E-2 

-1.65 

2 

22038-1-6437 

-5.52 

4.26E-7 

4.91E-5 

1.40E-4 

-0.46 

2 

22288-0001 

-4.03 

2.95E-4 

3.32E-4 

2.44E-2 

-1.87 

2 

22388-1-4419 

-1.86 

1.94E-2 

1.94E-2 

0.77 

-1.60 

2 

23078-1-7523 

-3.98 

8.76E-6 

1.18E-5 

2.08E-3 

-2.25 

2 

23393-1-4543 

-1.77 

5.05E-2 

5.08E-2 

0.72 

-1.15 

2 

23393-1-4543 

-1.86 

5.53E-2 

5.53E-2 

0.45 

-0.91 

2 


The formulas for calculating l 2 (t?)di? by elliptic functions can refer to Byrd & Friedman (1971) 
Similar studies which used elliptic functions can refer to |Kozai (1962)]~ Soderhjelm (1982) and 


Solovaya (2003)[ The remainder term Ti(t?) — T 2 (t^) is generally small and sometimes can be ig¬ 


nored. If Ti(t?) — T 2 ('d) can be ignored, 6 can be calculated analytically by elliptic functions. But 


here \Zi (z?) — I 2 ('d)] d'd is considered by simple Newton-Cotes integration formula to make a bet¬ 
ter approximation. 0 can be solved approximately by an iterative method. The three angular variables 
■^ 1 , 90 , hi can be integrated also by simple Newton-Cotes integration formula simultaneously. Another 
simplification is that the implicit Zeipel transformations from the averaged variables to the osculating 
elements can be turned into explicit. We call this model as MC-model. 


We now apply this model to 25 real triple stars with determined dynamical state (component masses 
and kinematic parameters). The results are listed in Table[T]including system name, order of magnitude 
of the perturbation (log^g )^ the RMSEs of M-model, K-model and MC-model, the ratio of 

the RMSE of MC-model to that of K-model (log^g ) and the type of motion. According to this table, 

the accuracy between M-model and MC-model is comparable. Eor all these stars, the RMSE of MC- 
model in comparison with the K-model’s, is reduced significantly. Indeed, for ^ 60% stars, the RMSEs 
are reduced by more than one order of magnitude. To show more details, we take WDS 02022H-3643 
as an example. Erom the N-model, the deviations of component positions calculated by M-model, MC- 
model and K-model, respectively, are shown in Eig. |2] Erom this figure, we know that the performance 
of MC-model is almost as good as M-model’s. When compared with K-model, the model accuracy is 
significantly improved and the applicable time span is significantly increased. 

As we all know, one of the important factors decide the quality of dynamical state determination is 
the accuracy of the dynamical model. In order to show the improvement in this respect brought by the 
high accuracy MC-model, we apply both this model and K-model to two systems, WDS 20396H-0458 
(HIP 101955, type 1) and WDS 00325H-6714 (HIP 2552, type 2). 
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Fig. 2 From the N-model, the deviations of component positions of WDS 02022H-3643 cal¬ 
culated by using M-model, MC-model and K-model, respectively. 


Two kinds of observations, relative position data (RPD) and the Hipparcos Intermediate Astrometric 
Data (HIAD) are used in the fitting. RPD are extracted from the Washington Double Star (WDS) 
Catalog (IMason et al. 200 11) . and the Fourth Catalog of Interferometric Measurements of Binary Stars 
( Hartkopf et al. 200T| |. HIAD are the abscissa residuals with respect to a reference point, the abscissa of 
which is calculated from a given solution. HIAD are read from the resrec folder on the catalogue DVD of 
|Leeuwen(2007)| With these observational data, the maximum likelihood estimate of model parameters 
is obtained by minimizing the objective function (x^) 




))= 


( 8 ) 


where yi is the observational quantity, y{xi] ui • • • om) is the corresponding calculated value according 
to the model parameters oi • • • om- We use the Bounded Variable Least Squares (BVLS) algorithm 
dLawson & Hanson 19951 ) to minimize the x^- 

HIP 101955 is a nearby low-mass triple star. There are 15 RPD points spanning from 1998 to 2008 
of inner orbit, 46 points from 1934 to 2008 of the outer one, and 91 HIAD in reference to a solution 
with 5 parameters. In the previous determinations of the dynamical state, the Kepler’s two-body motion 
model is applied separately to the inner {Aa, Ah} and the outer {Am, B} where Am is the center-of- 
mass of the inner binary AaAb ( Malogolovets et al. 2007|l. The resul ts are collected in the Sixth Catalog 
of Orbits of Visual Binary Stars(ORB6) ( [Hartkopf & Mason 2014| l, where the inner and outer orbits 
are roughly evaluated as good and reliable, respectively, according to the orbital coverage of the ob¬ 
servations. Because more observations are added, we firstly also use the K-model to fit observations. In 
comparison with the previous results, the is found to be reduced by ~ 66%. When the fitting model is 
replaced by MC-model, the x^ is further reduced by ^ 44%. Therefore, we conclude that using high ac¬ 
curacy MC-model, the fitting result is significantly better than the previous K-model’s results. Using the 
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Fig. 3 The fitting result of HIP 101955. 



X(arcsec) 



Fig. 4 The fitting result of HIP 2552. The open circle is a discarded point. 


fitted dynamical state parameters, the RMSEs of MC-model and K-model are calculated during the for¬ 
ward 100 years, that is, from 2008 to 2108. The RMSE of MC-model in comparison with the K-model’s, 
is significantly reduced by more than 80%, from 35.9mas (K-model) to ~ G.Omas (MC-model). This 
result shows that though starting with the same initial condition, for HIP 101955, the K-model can not 
be used to predict the component positions. 

Eor HIP 2552, there are 16 RPD points spanning from 1989 to 2005 of inner orbit, 75 points from 
1923 to 2010 of the outer one, and 151 HIAD in reference to an acceleration solution with 7 parameters. 
The inner and outer orbits were provided by [Docobo et al. (2008)| and are evaluated as good and inde¬ 
terminate by ORB6. K-model is also firstly used to fit the observations. In comparison with the previous 
fitting results, the is reduced by ~ 42%. When the fitting model is replaced by MC-model, though 
the is not significantly reduced, the RMSE is reduced from 10.5mas which is calculated by K-model 
to 0.74mas by MC-model. Using the fitted dynamical parameters, during the forward 100 years, the 
RMSE of K-model is 29.8mas while ^ S.Omas of MC-model. Therefore, K-model is also not suitable 
to predict the component positions for HIP 2552. 

We plot the fitted trajectories of HIP 101955 and HIP 2552, respectively, in Pig. |3] and Pig. |4] In 
these two figures, the filled circles are the RPD used in fitting, solid curves represent the previous double 
two-body model while the dotted curves are the fitted trajectories calculated using the MC-model. The 
trigonometric curves represent the N-model. As shown in the two figures, the difference between MC- 
model and the N-model is small enough to be ignored. The fitted dynamical state parameters and their 
Icr errors are listed in Table|2]and[3 
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Table 2 The fitted dynamical masses and kinematic parameters of HIP 101955. 


parameter 




unit 

M 

0.786±0.11 

0.493±0.11 

0.516ih0.21 

Mq 

^Ab 

-0.0598±0.0050 

0.127±0.0050 

-0.0188±0.022 

arcsec 

'f'B 

-0.188±0.0040 

0.173±0.0038 

0.902ih0.025 

arc sec 

'^Ah 

-0.102±0.0062 

-0.206±0.016 

0.1Ilih0.039 

arcsec/yr 


0.0367±0.0027 

-0.174±0.0066 

0.0669±0.016 

arcsec/yr 


Table 3 The fitted dynamical masses and kinematic parameters of HIP 2552. 


parameter 




unit 

M 

0.389±0.038 

0.0969±0.038 

0.177 ±0.212 

M0 

'^Ab 

-0.0614±0.047 

-0.298±0.029 

0.290±0.032 

arcsec 

VB 

-4.029±0.016 

0.609±0.015 

-0.318±1.8 

arcsec 

^Ab 

0.235±0.015 

-0.0331±0.013 

0.000668±0.025 

arcsec/yr 

VB 

0.0478±0.0059 

-0.0556±0.0029 

0.0455±0.0097 

arcsec/yr 


5 CONCLUSION AND DISCUSSION 

Marchal’s first-order analytical solution is implemented and a more efficient simplified version is applied 
to real hierarchical triple stars. The results show that the proposed first-order model is preferable to 
the classical double two-body model both in fitting observational data and in predicting component 
positions. 

As pointed out in section[3 there are a few cases to which the M-model doesn’t apply, because of the 
inadequacy of the Delaunay elements. For these cases, Poincare elements should be used instead. There 
are also a few cases when the first-order perturbations are very small in the time span of observations, 
but its maximum value over the whole period of the outer orbit is too large to apply M-model. For 
these cases, our preliminary studies show that it is possible to give a suitable first-order solution without 
resorting to averaging over the outer orbit. 
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